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Abstract — In the context of the compressed sensing problem 
we propose a new ensemble of sparse random matrices which 
allow (i) to acquire/compress a po -sparse signal of length TV in a 
time linear in TV and (ii) to perfectly recover the original signal, 
compressed at a rate a, by using a message passing algorithm 
(Expectation Maximization Belief Propagation) that runs in a 
time linear in TV. In the large TV limit, the scheme proposed here 
closely approach the theoretical bound po = a, and so it is both 
optimal and efficient (linear time complexity). More generally, 
we show that several ensembles of dense random matrices can 
be converted into ensembles of sparse random matrices, having 
the same thresholds, but much lower computational complexity. 

I. Introduction 

Compressed sensing is a framework that enables to recover 
an TV-dimensional sparse signal s = (si) from M(< TV) linear 
measurements of its elements, y = Fs, exploiting the prior 
knowledge that s contains many zero elements 0]. A simple 
consideration guarantees that the ^o- recover y 



s = argmin ||a;||o subj. to y = Fx, 



(1) 



where ||a?||o denotes the number of non-zero elements in x, 
is theoretically optimal in terms of minimizing the number 
of measurements M necessary for perfectly recovering any 
original signal s. However, carrying out the £o- recover y f° r 
a general measurement matrix F is NP-hard. For avoiding 
such computational difficulty, an alternative approach, the t\- 
recovery 



s = argmin 1 1 x \ | i subj. to y = Fx, 



(2) 



where ||a;||i = X)j=i l x il> ^ s widely employed, as (f2]l is 
generally converted to a linear programming problem and 
therefore signal recovery is mathematically guaranteed in an 
0( TV 3 ) computational time by use of the interior point method. 
Nevertheless, the 0(TV 3 ) cost of computation can still be 
unacceptable in many practical situations, and much effort is 
being made for exploring more computationally feasible and 
accurate recovery schemes lO-Q. 

Among such efforts, a recovery scheme recently proposed 
by Krzakala et al. Q is worth of attention. Their scheme 
basically follows the Bayesian approach. Namely, the signal 
recovery problem is formulated as that of statistical inference 
from the posterior distribution 

6{Fx - y)P{x) 



P(x\F,y) 



Z(F,y) 



(3) 



where Z(F,y) is a normalization factor imposing the con- 
dition J dxP(x\F,y) = 1 and a component-wise prior 

distribution P(x) = YiiLiii^ — p)^( x i) + P4>{ x i)\ i s as ~ 
sumed. p and <fi(x) represent the density of non-zero signal 
elements and a Gaussian distribution, respectively. Exactly 
inferring s from ([3]) is NP-hard similarly to ([1). However, 
by employing the belief propagation (BP) in conjunction with 
the expectation-maximization (EM) algorithm for estimating p 
and the parameters of 4>(x), they developed an approximation 
algorithm, termed EM-BP, which overcomes the recovery 
performance of (O with an only 0(N 2 ) computational cost. 
Furthermore, they showed that, by employing a peculiar type 
of "seeded" matrix F, the threshold value of the compression 
rate a — M/N of EM-BP, above which the original signal 
is typically recovered successfully, can approach very close 
to that of the ^-recovery, (x s -embp = Po, where po is the 
actual signal density of s and s — EM BP stands for 'seeded 
EM-BP'. The seeded matrix is composed of blocks along 
the diagonal densely filled with Gaussian random variables. 
It is important to remember that this result is achieved for 
the first time with an approach different from the £o- recover y, 
being the threshold for the l\ -recovery much higher than the 
optimal one: ai 1 > po- As the optimality of the £q -recovery is 
guaranteed for the EM-BP, this indicates that this scheme can 
practically achieve the theoretically optimal threshold of the 
signal recovery with an 0(N 2 ) computational cost. However, 
the optimality of their scheme in terms of the computational 
complexity is still unclear; there might be a certain design of 
the measurement matrix F that makes it possible to further 
reduce the necessary computational cost keeping the same 
signal recovery threshold. 

The purpose of the present work is to explore such a 
possibility. For this, we focus on a class of matrices that are 
characterized by the following properties: 

• sparsity: The matrix F has only O(l) non-zero elements 
per row and column. This implies that the measurements 
can be performed in a time linear in the signal length. 
This is highly preferred for practicality, given that such an 
operation typically needs to be done in real time, during 
data acquisition. 

• integer values: The matrix elements are not real valued, 
but take small integer values. This means that an opti- 
mized code for the measurements can work with bitwise 



operations, thus achieving much better performance with- 
out any loss of precision. 
• no block structure: The block structure used in [7 1 may 
not be necessary for reaching the optimal threshold. As 
an alternative possibility, we study a structure made of 
a square matrix in the upper left corner (the seed) plus 
a stripe along the diagonal. This structure is much more 
amenable for analytic computations, since it corresponds 
to a one-dimensional model homogeneous in space. 
The use of sparse matrices for compressed sensing has al- 
ready been proposed in several earlier studies ll8l- lfl2l . Among 
them, our approach is particularly similar to that of lfl2l in the 
sense that the both sides are based on the Bayesian frame- 
work and use integer-valued sparse measurement matrices. 
Nevertheless, the two approaches differ considerably in the 
following two points. Firstly, we adopt EM-BP, which updates 
only a few variables per node for the signal recovery, while 
the recovery algorithm of [12] involves functional updates 
and needs significantly more computational time than ours. 
Secondly, we make much effort to explore a simple design of 
F that achieves a nearly optimal recovery performance while 
the problem of the matrix design is not fully argued in [ 12 1. By 
carrying out extensive numerical experiments in conjunction 
with an analysis based on density evolution [13], we show 
that a threshold close to the theoretical limit a — po can be 
achieved by the matrices with the above properties with almost 
linear computational cost in both measurement and recovery 
stages. 

This paper is organized as follows. In Section [TT] the EM- 
BP algorithm is briefly explained and the results for dense 
matrices are summarized. The algorithm is applied to homoge- 
neously sparse matrices in Section [Til] and to structured sparse 
block matrices in Section [IV] A new type of "striped" sparse 
matrix without blocks is introduced in Section [V] The last 
Section summarizes the work focusing on its importance for 
practical use and on perspectives. 

II. Expectation Maximization Belief Propagation 

The new algorithm based on BP in conjunction with EM 
proposed in Ref. |7| starts from Eq. OJ. In order to solve it 
with BP, 0(MN) messages for the probability distributions 
of the variables Xi are constructed in the following way: 

1 



derives from an expansion in small Fai, and in the dense 



(Xi) 



Z' 1 



m^^Xi) = — — [(1 - p)5(xi) + p</>(xi)^ Y[ 
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where Z 1 ^^ and Z^ % are normalization factors. This EM- 
BP equations are very complicated because the messages 
are distribution functions. In order to make them simpler, 
the messages can be approximated assuming that they are 
Gaussian, thus obtaining the equations for the mean 
and the variance Wj-j.^ of mi^^Xi). This approximation 
was firstly introduced for sparse matrices in Ref. [14| and 
it becomes exact asymptotically if F is dense. In fact it 



case F^i = 0(1/ yN). Supposing that the elements of the 
original signal follow a Bernoulli-Gaussian distribution with 
parameters po, xq and ao, the update rules for the messages 
are the following: 

a i— = fa I ^ ' A^^ S ' B^^. t 
\7#M 7#M 
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(4) 



where f a and f c are some analytical functions depending on 
the parameters p, x and a. For details, see Ref. Q. 

In general, the original p , xq and ao are not known, but 
one can use EM to derive the update rules for them, using the 
property that the partition function 

Z(p,x,a)= / dxP(x)6(y - Fx) 



is the likelihood of the parameters (p, x, a) and is maximized 
by the true parameters po, xq and ao. Thus, after the update 
of all the messages, the inferred parameters of the original 
distribution are updated following these rules: 



x <- 
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and V; 



_ 1 - e 2(l/cr-' + U i ) 2a 



with Ui = 2 7 A ,« - , , _ . 

If the algorithm converges to the correct solution, a; = s, 
and Vi = 0. 

To reduce the number of messages from 0(NM) to 0(N), 
one can notice that in the large N limit, the messages and 
are nearly independent of p. Thus equations involving 
only a variable per each measurement node and a variable per 
each signal node can be derived, being careful to keep the 
correcting Onsager reaction term as for the TAP equations 
in statistical physics H151 . This method was introduced in 
the context of compressed sensing in Ref. [6| and is called 
approximated message passing (AMP). 

In general, the correct distribution of the original signal 
is unknown. However, in Ref. Q it is demonstrated that if 




Fig. 1. Top: Probability of perfect recovery versus the signal sparsity po 
using sparse matrices with a = 0.5 and K = 20. The threshold is the same 
as with dense matrices. Bottom: Probability of perfect recovery computed 
with density evolution has the same threshold. Here TV is the population size. 



a > po, the most probable configuration of x with respect to 

P{ x ) = rii=i[(l — P)3{ x i) + P$( x i)] w i tn P < 1> restricted 
to the subspace y = Fx, is the original signal s, even if the 
signal is not distributed according to P(x). So our choice of 
a Gaussian distribution for (j>{x) should be perfectly fine also 
even if the original signal has a different distribution. 

If a dense matrix is used, the free entropy $(-D) at fixed 
mean square error D = (1/N) — s *) 2 can he 

computed. For a > po, the global maximum of the function 
<£>(£)) is at D = 0, that corresponds to the correct solution. 
However, below a certain threshold a < a bp that depends on 
the distribution P(s), the free entropy develops a secondary, 
local maximum at D ^ 0. As a consequence, the EM- 
BP algorithm can not converge to the correct solution for 
po < a < asp, because a dynamical transition occurs. 
Nonetheless, the threshold olb p is lower than ai x . 

III. EM-BP WITH A SPARSE MATRIX 

First of all we want to verify if the use of a sparse matrix 
can reach the same results than the use of a dense one. For 
a sparse random matrix the AMP equations can not be used, 
thus we use the update rules in Eq. dU, for the inference of the 
original signal s. In particular, we choose the matrix F to have 
only K — 0(1) elements different from zero in each row and 
H = aK = 0(1) elements in each column, extracting them 




Fig. 2. Stability check for the correct solution under the EM-BP message 
passing algorithm. Starting the recovery process with a sparse matrix from an 
initial condition differing less than A from the correct solution, the latter is 
recovered as long as po < c*stab(A). In the limit A — > the stability limit 
Cstab(A) tends to the theoretical bound a (which is 0.5). 

from the distribution 

P(Ffii) = \S{F^ - J) + \S{F^ + J) (5) 

with J = 1. The use of the messages a^^x) and Vi- > . fl (x) 
instead of the AMP equations does not involve an extra cost in 
memory, because, due to the sparsity of the matrix, the number 
of the messages is O(N). In principle the messages mj^fa;) 
are not Gaussian if the matrix is sparse, so the use of only the 
two parameters a^^x) and Vi-y^x) is not exact. However, 
the convolution of K messages (with K = 20 in a typical 
matrix we use) is not far from a Gaussian and indeed we can 
verify a posteriori that this approximation is valid, because it 
gives good results. 

In all our numerical simulations we use a Bernoulli- 
Gaussian distributed signal and a compression rate a = 0.5. 

In Fig. [T] (top) the probability of a perfect recovery as a 
function of the sparsity of the signal po for different sizes is 
shown, applying the EM-BP algorithm using a sparse matrix 
with K — 20. The threshold for the perfect recovery in the 
thermodynamic limit (TV — > oo) is psp — 0.315, that is 
the same obtained in Ref. Q with a dense matrix. We can 
not analytically compute the free entropy Q(D) as in Q, 
because we are using sparse matrices and methods such as the 
saddle point one can not be used. However we can perform 
a numerical density evolution analysis, as shown in Fig. Q] 
(bottom) and the threshold is the same as the one computed 
with the matrices F. 

Next, we have verified that the correct solution is always 
the global maximum of < &(£ ) ) and it is locally stable up to 
a ~ po using EM-BP with a sparse matrix. Since we can 
not compute analytically the free entropy, we must resort to a 
numerical method. We start EM-BP with an initial condition 
very close to the correct solution: a°_^ M = s, + Si^, with 
Si-tn a random number uniformly distributed in [—A, A]. In 
this way we have verified (see Fig. |2]i that if A is sufficiently 
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Fig. 3. Pc(N) with a sparse, structured matrix with blocks, for different 
sizes N and values of L (see text). The thresholds for the l\ -recovery and 
for EM-BP without any structure are also drawn. For comparison with data 
in Fig. [3] we have used the same exponent a ~ 0.18 best fitting those data. 

small, the correct solution can be found up to a ~ po, as in 
the case of a dense matrix. 

For the algorithms based on the l\ minimization, it is known 
that the threshold with a sparse matrix is lower than that with a 
dense one. However, those algorithms are not optimal, because 
the correct solution disappears below the threshold ai t . In this 
sense, the EM-BP algorithm is optimal, because the global 
maximum of the free entropy is always on the correct solution. 
Thus one can expect that, if the rank of the sparse matrix is 
the same as that of the dense one, a similar threshold can be 
reached, as we have demonstrated numerically. 

Concluding this first section, we can say that the EM- 
BP algorithm of Ref. Q seems to reach the same threshold 
olbp, either using a dense Gaussian matrix or a sparse binary 
one. However the use of a sparse matrix is computationally 
much faster than a dense one. Moreover, the use of binary 
elements, instead of Gaussian real values, allows for a better 
code optimization and eventually for hard-wiring encoding of 
the compression process. 

IV. Block-structured sparse matrices 

To overcome the secondary maximum of the free entropy, 
in Ref. [7| the authors use a kind of structured block matrix 
that helps to nucleate the correct solution. The idea is that the 
correct solution can be found for the first variables, and then it 
will propagate to all the signal. This idea is similar to the so- 
called spatial coupling, that is very useful for many different 
problems [16|. With this trick the authors of Ref. [7| reach a 
perfect recovery for almost any a > po in the large N limit. 
Here we try to use a matrix with the same block structure, 
but sparsely filled. We divide the N variables into L groups 
of size N/L and the M measurements into L groups of size 
M p = ctpN/L in such a way that M = Yl P =i M p — aN and 

1 / L J2p=i a P = a - 

In this way the matrix F is divided in L 2 blocks, labeled 
with indices (jp, q). Each block is a sparse binary matrix with 



k elements different from zero for each row and h p = a p k 
elements for each column, distributed according to Eq. ||5), 
with J = J p>q . As in Ref. Q we choose J P}P -i — Ji, 
Jp.j> — 1, Jp,p+i = Ji and Jp. q = otherwise. The important 
ingredient to nucleate the correct solution is that in the first 
block a.\ = (M\/N)L > asp holds. For simplicity we 
choose ot\ = 1 and a p = (La — 1)/(L — 1) for p ^ 1. The 
dependence of the recovery success from the parameters J\ 
and J2 is strong and the best results for a = 0.5 are obtained 
around J\ = 4 and J2 = 1. We have chosen this two values 
for the following experiments also because we want to work 
with matrices whose elements are small integer values. 

Similarly to the dense case, the use of a sparse structured 
matrix with blocks allows to overcome the dynamical transi- 
tion at olbp and to nucleate the correct solution until a very 
close to pq. In Fig. [3] the mean critical threshold p c (N) is 
shown for different signal lengths at a fixed compression rate 
a = 0.5. On the x axis we have used the same scaling variable 
as in Fig. [5] and the best parameter a obtained from the fit of 
data in Fig. also interpolates quite well data in Fig. [3] In the 
thermodynamic limit, p c extrapolates to a value compatible 
with the optimal one, a, and it is certainly much higher than 
the thresholds for the l\ -recovery and for EM-BP without any 
structure. We have done also a density evolution analysis, that 
confirms this result. 

For each value of N and L, the mean critical threshold p c is 
computed as follows. We randomly generate a block structured 
matrix F with given N and L. We start with an original signal 
s sufficiently sparse, which is recovered by the algorithm; then 
we add non-zero entries to the signal, and check whether the 
new signal can be recovered by the algorithm; we go on adding 
non-zero elements to the signal, until a failure in a perfect 
recovery occurs. The previous last value for po is the critical 
threshold for the matrix F. The mean critical threshold is 
obtained by averaging over many different random matrices 
and signals, with the same values of N and L. The number 
of such random extractions goes from 10 3 for the largest N 
value up to 10 4 for the smallest N value. 

The values of (N, L, k) used for the simulations shown 
in Fig. E] are the following: (2250, 10, 9), (19000, 20, 19), 
(31200, 40, 39) and (49000, 50, 49). We need to increase both 
N and L if we want to obtain good results in the thermo- 
dynamic limit. However, if we change L, we must change k 
too. Indeed, in order to have the same number of elements per 
row and column for each of the L 2 blocks, we must satisfy 
the conditions: (N/L)/M p = k/h p with k and h p integer 
valued. We have used the smaller possible value for k, that is 
k = L — 1. The fact that it is impossible to keep k constant 
while increasing L, implies that these kind of block-structured 
matrices always become dense in the thermodynamic limit. 
This is a limitation of the block structure that we want to 
eliminate with the matrix proposed in the following Section. 

V. And without blocks? 

The matrix proposed in Ref. [7| is not the only one that 
allows to reach the optimal threshold. In Ref. ifTTl the use of 
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Fig. 4. Nearly one-dimensional sparse matrix with a first squared block of 
size L X L, and non-zero elements in the stripes around the diagonal, can 
achieve compression and perfect recovery close to the theoretical bound in 
lineal' time. 



Fig. 5. pc(N) with a sparse, striped matrix, as described in Section IV1 for 
different sizes (from N = 2000 to N = 40000). The thresholds for the 
^l-recovery and for EM-BP without structure are drawn for comparison. The 
best fitting parameter is a ~ 0.18 and leads to an extrapolation of p c (N) in 
the thermodynamic limit compatible with a. 



other good dense, block-structured matrices is analyzed. How- 
ever the block-structure is not so simple to handle if one wants 
to do analytical calculations in the continuum limit. Moreover 
in making these block-structured matrices sparse one has to be 
careful to find right values for L, M, N, M p , a p , k. For these 
reasons we want to understand if the block structure is crucial, 
and, if not, we want to eliminate it. 

We try a different structured sparse matrix (see Fig. 0J, that 
we call striped matrix. It has one sparse square block of side 
L on the top left of the matrix with K = 0(1) elements for 
each row and column extracted from <(5j with ,7=1. This is 
fundamental to nucleate the correct solution. Apart from this 
first block, the residual compression rate is a' — A ^Zl - Then 
we construct a one-dimensional structure, around the diagonal 
of the remaining matrix. For each column c > L, we put 2Ka' 
non-zero elements, again extracted from (0, randomly placed 
in the interval of width 2La' around the diagonal. One element 
with J = 1 is always placed on the diagonal (actually on the 
position closest to the diagonal). For the remaining elements 
we use the following rules. If the element is at distance d < 
La'/3 from the diagonal, we use J = 1. Otherwise, if its 
distance is d > La' /3, we use J = J\ below the diagonal 
and J — J2 above the diagonal. In this way the number of 
elements per column is constant, while the number of elements 
for row is a truncated Poisson random variable with mean 2K : 
indeed there are no empty rows thanks to the rule of placing 
the first element of each column closest to the diagonal. When 
constructing the matrix we apply exactly the same rule to each 
column, but on the last L columns it may happen that a non- 
zero element has row index larger than M: these elements 
are then moved below the first squared matrix by changing 
row and column indices as follows: r <— r — (M — L) and 
c <- c — (N — L). 

In this way we have some kind of continuous one- 
dimensional version of the block-structured matrix discussed 
in the previous Section. Within this striped matrix ensemble 



the thermodynamic limit at fixed matrix sparsity can be taken 
without any problem, by sending N, L —> 00 at fixed L/N 
and fixed K = 0(1). In Fig. [5] we show the mean critical 
threshold reached by using striped matrices with a fixed ratio 
L/N = 1/50 (the same used in the plot of Fig. |4l) and different 
signal lengths. Perfect decoding up to p c is again achieved by 
using the EM-BP algorithm. We have extrapolated the p c (N) 
data to the thermodynamic limit by assuming the following 
behavior in the large N limit 

p c (N) = p c (oo) - bN- a (6) 

Data in Fig.|5]have been plotted with the best fitting parameter 
a ~ 0.18 and the extrapolated value p c (oo) is perfectly 
compatible with the theoretical bound a. 

So we can conclude that the important ingredients to reach 
optimality is not the block structure, but the nearly one- 
dimensional structure, associated with an initial block with 
ai > asp to nucleate the correct solution. 

It is worth noticing that the corresponding statistical me- 
chanics model for these striped random matrices is a one- 
dimensional disordered model with an interaction range grow- 
ing with the signal length as in a Kac construction. Models of 
this kind are analytically solvable, and show very interesting 
results |[T8ll . 

The use of our striped sparse matrices allows for a great 
reduction in the computational complexity. Indeed, both times 
for measurement and recovery grow linearly with the size 
of the signal if sparse matrices are used, while they grow 
quadratically if dense matrices are used. In Fig.|6]the times for 
the measurement/recovery of a signal are shown for different 
signal lengths N. For this test, we used dense block-structured 
matrices and sparse striped matrices. The number of EM-BP 
iterations to reach the solution is roughly constant for different 
N. A quadratic fit for the dense case and a linear fit for the 
sparse one perfectly interpolate the data. 
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Fig. 6. Actual time (in seconds) for the recovery of a signal with dense and 
sparse matrices for different data length N. The data are fitted respectively 
by a quadratic and a linear function. 

VI. Conclusions and future developments 

We have introduced an ensemble of sparse random matrices 
F that, thanks to their particular structure (see Figure^, allow 
to perform in a linear time the following operations: 

(i) measurement of a po-spaise vector s of length N by a 
linear transformation, y = Fs, in a vector y of length otN; 

(ii) perfect recovery of the original vector s by a message 
passing algorithm (Expectation Maximization Belief Propaga- 
tion) for almost any parameter satisfying the theoretical bound 
Pa < a. 

The main ingredients of these striped sparse matrices that 
allow to reach such good performances are the presence of 
a 'seeding' sparse square matrix in the upper left corner that 
nucleates a seed for the right solution and the one-dimensional 
structure along the diagonal that propagates such an initial seed 
to the complete right solution. Both the seeding and the one- 
dimensional structure have been already used in the past Q, 
JT6), but in our new ensemble the matrices are sparse and this 
permits to perform all the operation in a time linear in the 
signal length. 

We have also checked that sparse matrices perform as well 
as dense ones in case of block-structured matrices and for 
matrices with no structure at all. 

Apart from the compressed sensing case, the sparsity of the 
matrix or equivalently the linear time complexity is essential 
in several applications J5) 

In data streaming computing, one is typically interested in 
doing very quick measurements in constant time. For example, 
willing to measure the number of packets Sj with destination 
i passing through a network router, it is not possible to keep 
the vector s because is generally too long. Instead, a much 
shorter sketch of it, y = Fs, is measured, such that the very 
sparse vector s can be recovered from y. The sparsity of the 
matrix F is essential in order to be able to update the sketch 
y in constant time for each new packet passing through the 
router. 

Another interesting application is the problem of group 



testing, where a very sparse vector s 6 {0, 1}^ is given and 
one is interested in performing the fewest linear measurements, 
y = Fs, that allow for the detection of the defective elements 
(Si = 1). In this case the sparsity of the matrix F is required 
by the experimental constraints: only if the tested compound 
= Y^,i F^Si, is made of very few elements of s the linear 
response holds and non-linear effects can be ignored. 

However in the more general case, one does not observe 
directly the sparse signal s but rather a linear transformation 
of it, x = Ds, made with a dictionary matrix D (which 
is typically a Fourier or wavelet transformation, and thus 
is a dense matrix). In this more difficult case, one would 
like to design a sparse measurement matrix A such that the 
measurement/compression operation, y = Ax, is fast and the 
resulting observed data y is short thanks to the sparseness 
of s. The conflicting requirement is to have a fast recovery 
scheme, because now to recover the original signal one should 
solve s = argmin||s||o subject to y = (AD)s, where 
AD is typically dense (e.g. in case of Fourier and wavelet 
transformations). So a very interesting future development of 
the present approach is to extend it to this more complex case. 
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